Kinetic Temperatures for a Granular Mixture 



Steven R. Dahi and Christine M. Hrenya 
Department of Chemical Engineering, University of Colorado, 
Boulder, CO 80309 

Vicente Garzö 

Departamento de Fisica, Universidad de Extremadura, E- 06071 

Badajoz, Spain 

James W. Dufty 
Department of Physics, University of Florida, Gainesville, FL 32611 

(February 1, 2008) 

Abstract 

An isolated mixture of smooth, inelastic hard spheres supports a homoge- 
neous cooling state with different kinetic temperatures for each species. This 
phenomenon is explored here by molecular dynamics simulation of a two com- 
ponent fluid, with comparison to predictions of the Enskog kinetic theory. 
The ratio of kinetic temperatures is studied for two values of the restitution 
coefficient, a = 0.95 and 0.80, as a function of mass ratio, size ratio, compo- 
sition, and density. Good agreement between theory and simulation is found 
for the lower densities and higher restitution coefficient; significant disagree- 
ment is observed otherwise. The phenomenon of different temperatures is also 
discussed for driven systems, as occurs in recent experiments. Differences be- 
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tween the freely cooling state and driven steady states are illustrated. 
PACS number(s): 45.70.Mg, 05.20.Dd, 51.10.+y, 47.50.+d 
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I. INTRODUCTION 



The dissipative nature of granular media is captured by an idealized fluid of smooth, 
inelastic hard spheres. When isolated and homogeneous such a system rapidly approaches a 
homogeneous cooling state (HCS) for which ali time dependence of the distribution function 
occurs through the temperature. The latter, defined in the usual way via the average kinetic 
energy, decays in time ("cooling") due to the inelastic collisions. The existence of the HCS 
and associated cooling rate is well established for a one component system by theory [[T], 
Monte Carlo simulation 0, and molecular dynamics simulation ||. Recently, it has been 
shown from the Enskog kinetic theory that a mbcture of inelastic hard spheres also has a HCS 
under the same conditions Q. The condition that ali time dependence occurs through the 
temperature requires that the cooling rates for the kinetic temperatures for each species must 
be the same. It follows directly that the kinetic temperatures are different for mechanically 
different species, reflecting a violation of the equipartition theorem valid for elastic collisions. 
A prediction for the ratio of temperatures in a binary mbcture as a function of mass ratio, size 
ratio, composition, density, and restitution coefficients was obtained from an approximate 
solution to the Enskog equations. The accuracy of this approximate result has been recently 
confirmed by Monte Carlo simulation of the Enskog equations ||. 

The objective here is to demonstrate the phenomenon of a HCS and two temperatures in 
a broader context by molecular dynamics (MD) simulation for a binary mixture of inelastic 
hard spheres. MD simulation avoids any assumptions inherent in the kinetic theory or 
approximations made in solving the kinetic equations. It is shown here that MD simulation 
supports the existence of a HCS for mixtures with different kinetic temperatures for each 
species but a common cooling rate. The dependence of the temperature ratio on mechanical 
properties and state conditions is found to be in good agreement with predictions of the 
Enskog kinetic theory except at high density and strong dissipation. In the latter case, 
significant quantitative deviations from the Enskog theory are observed but the concept of 
a HCS and two temperatures is preserved. 
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The HCS can be given a time independent representation by transformation to suitable 
dimensionless variables |7|||. In this form it is similar to the steady state obtained for 
homogeneously driven granular fluids. The latter are obtained by adding stochastic sources 
to the kinetic equation or MD dynamics to do work on the system that compensates for 
collisional cooling. The resulting homogeneous steady state is qualitatively the same as 
the dimensionless HCS, but the quantitative differences are expected to make it closer to 
locally driven steady states observed in experiments on vibrated granular media. Studies 
of driven states have been extended to mbrtures both theoretically |§ and experimentally 
n^Ijfl . Comparisons of the temperature ratio for the HCS mbcture and that for two types of 



homogeneously driven mbctures is given below. Their relationship to a locally driven system 
is discussed also. 

The plan of the paper is as follows. In Sec. |H], we show that the Liouville equation for a 
binary granular mbrture supports a scaling solution describing the HCS. A transformation 
to dimensionless variables allows to get the (constant) temperature ratio 7 = T 1 (t)/T 2 (t) 
in terms of the parameters of the mbcture. An approximate evaluation of the temperature 
ratio can be made from the Enskog kinetic theory, as is shown in Sec. fTT[ In Sec. |TV], 
the Enskog predictions are compared with those obtained from MD simulations. Such a 
comparison shows a quite good agreement for the lower densities considered but significant 
discrepancies between theory and simulation appear for high density and strong dissipation. 



The existence of two temperatures in driven granular mixtures []9| jT^ , p~3[] and its possible 
connection with recent experiments is analyzed in Sec. [V|. Finally, the paper ends in Sec. 



VI with a brief discussion on the relevance of the results presented here. 



II. HOMOGENEOUS COOLİNG STATE FOR A MIXTURE 

The system considered is a binary mixture of N\ and N2 smooth hard spheres of masses 
mı and 1712, and diameters <Tı and o~2- In general, collisions among ali pairs are inelastic 
and are characterized by three constant restitution coefficients a^-, which can be different 
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for the three types of pair collisions. The state of the system at time t is specified by the 
N = Nı + N 2 particle phase space density p(T, t) which is a solution to the Liouville equation 



I4]| . In ali of the following, attention is restricted to spatially homogeneous states. In this 
Section it is further assumed that the system is isolated. The properties of primary interest 
are the overall temperature T (t) associated with the total kinetic energy, and the partial 
temperatures T, (t) associated with the kinetic energies of each species. They are defined as 

2 3 / Ni 1 \ 

T(t) = J2x l T l (t), Nl T l {t) = (Y.-m l vl-t). (1) 

i=l \/x=l / 

The brackets denote a phase space average över the state of the system at time t and 
Xi = Ni/N is the composition. The time dependence of T (t) and T, (t) follows from the 
Liouville equation which gives P JT4]]1"5|| 



T-%T = -(, T- l d t Ti = -Q, (2) 

where Q is the cooling rate associated with the partial temperature Tj and ( is the total 
cooling rate. They are given by 

1 2 

C = 7p^2xiTiCi, (3) 

1 i=l 



m 



J dv&l J dv 2 J d,T 12 Tij (r 12 , vı, v 2 ) f-f (r 12) v 1; v 2 ,t) . (4) 



"in {T i ._j 

Here, rti is the number density of species i, rı 2 is the relative position of the two particles, 

(2) 

and flj (ri2, vı, v 2 , t) are the reduced two particle distribution functions for a particle of 
type i and one of type j, obtained from p(T, t) by integrating över degrees of freedom for ali 



other particles. The binary collision operators are defined by E, 15 



Tij (vı, v 2 , r 12 ) = crj / da Q{â ■ gı 2 )(o- ■ gı 2 ) a y 2 (J(rı 2 - (r)b^ - S(r 12 + a) (5) 

where <7y = (<Tj + <jj) /2, â is a unit vector directed along the line of centers from the sphere 
of species i to that of species j at contact, 9 is the Heaviside step function, and g 12 = vı— v 2 . 
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Also, b^ 1 is a substituting operatör, 6 İJ 1 F(g 12 ) = F (b^g^J, which changes any function of 
V! and v 2 to the same function of the restituting velocities v' x and v 2 : 

v [ = Vl - fiji (l + a^ 1 ) (d- • gıa)^, v 2 = v 2 + //y (l + c^ 1 ) (5- • g 12 )a- (6) 

where //jj = mı/ (m.j + m^). Upon writing Eqs. (f|) and (|) we have taken into account that 
for an homogeneous system the spatial dependence of fy occurs only through r 12 . 

In general, ali three temperatures and associated cooling rates will be different and 
depend on the initial preparation. The time evolution of the ratio of the two partial tem- 
peratures 7 (t) = Tı(t)/T 2 (t) follows from the second equality of Eq. (|2]): 

ö t ln 7 = C 2 -Cı- (7) 

For a system with elastic collisions, p(T, t) rapidly approaches to the Gibbs distribution with 
a single constant temperature. This requires Cı = C2, and Tj cx T in the Gibbs state. The 
form of the velocity distribution functions and the constancy of the temperature then gives 
Tı = T 2 = T and Cı = C2 = 0. This equality of the temperatures is the equipartition theorem 
for classical statistical mechanics. The vanishing of the cooling rates is a consequence of the 
system approaching towards a steady state. 

If the collisions are inelastic, the system stili approaches rapidly a special state known 
as the homogeneous cooling state (HCS). As with the Gibbs state, the velocities scale with 
the temperature for a dimensionless universal distribution of the form 

Pkcs(T,t) = Mt)}- 3N pl cs ({v* J ,vn). (8) 

Here r* 3 - = r^/İ denotes the dimensionless relative coordinate for particles i and j, and 
i is some appropriate characteristic length scale such as the mean free path. The di- 
mensionless velocities v* = Vj/ı>o(t) are scaled relative to the thermal velocity defined by 
vo(t) = y / 2T(t) (mı + m 2 ) Jm^m^. This scaling has the same consequences as described 
above for elastic collisions: Ti(t) oc T (t), j(t) — > constant, and so Cı = C2- However, Eq. 
(§) is not a steady state since the cooling rates Q do not vanish. Also, the form of p* hcs 
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is not the same as for the Gibbs state so there is no a priori reason to expect that the 
temperatures should be equal ||. In fact, as indicated below, they are equal only in the 
limit of mechanically equivalent particles or elastic collisions. 

The simplest test of the evolution towards a HCS with the assumed velocity scaling is 
the approach of 7 (t) to a constant value. This is illustrated in Fig. [I] from MD simulation 
using 1000 particles of the same size and composition, but with the mass ratio mı/777.2 = 8. 
The total solid volume fraction is = 0.1 and ali coefficients of restitution are equal. Here, 
= 0ı + 02, where 

0i = ^riicrf (9) 

is the species volume fraction of the component i. We consider two values of the resti- 
tution coefficient: a = 0.8 and a — 1. We observe that in the elastic case (a = 1) the 
mixture approaches the Gibbs state with the temperature ratio j(t) — > 1, as expected from 
equipartition. In the inelastic case (a = 0.8), 7(i) approaches a constant value (7^2 with 
fluctuations less than 5%) after about 10 collisions per partide. It is seen that the the HCS 
for inelastic collisions is approached on the same time scale as the Gibbs state for elastic 
collisions. Further details of the MD simulation are discussed below. 

III. ENSKOG KINETIC THEORY 

The kinetic temperatures defined by Eq. ([!]) can be given in an equivalent form in terms 
of the one partide reduced distribution function /j (v, t) as 

T,{t) = \- I dvv 2 f^(v,t). (10) 

The one partide reduced distribution function (v, t) obeys the exact first BBGKY hier- 
archy equations 
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To be more specific about the dependence of the temperatures on the parameters of the 
mbrture it is sufficient to specify the reduced distribution functions fj?' (rı 2 , Vı, v 2 , t) in Eq. 
(|TT|). This also determines the cooling rates from Eq. (^). These distribution functions occur 
only in the combination Tijfj?' , so knowledge of f^ is required only for pairs of particles 
at contact and only on the pre-collision hemisphere. A practical approximation for these 
conditions is obtained by neglecting velocity correlations and expressing the two particle 
distribution functions in terms of the single particle distribution functions 

(rı 2 , Vl) v 2) t) - /« ( Vl , t) /« (v 2 , t) Xij (r 12 , t) . (12) 

The single particle distributions are independent of position since only homogeneous states 
are considered here. The spatial correlation function Xij ( r i2,£) is evaluated at contact and 
its choice is given below. The one particle distribution functions can be determined by using 
this same approximation in the exact first BBGKY hierarchy equations, which becomes 



SjfVı,*) = E J ü N-^^J , (13) 

3=1 



2 



where Jij[fl , /i 1 ^] is the Enskog collision operatör 0]. These are now closed equations for 



fi and constitute the Enskog kinetic theory for the granular binary mbrture. 
For the HCS the scaling form (§) implies a similar scaling form for f ; (vı,t) 

ri 1 \v,t) = n i vû 3 (t)f*(v*). (14) 

Furthermore, \ij ( r i2 = <^ij,t) —> Xij = constant since ali time dependence occurs through 
the velocity scaling. For practical purposes, and to agree with the equilibrium limit for elastic 
collisions, Xij is taken to be the equilibrium pair correlation function. A good approximation 
is given by the Carnahan-Starling form [|Î6J 

v = -J— + 3 İ + 1 ? ( ^iY (i*) 

XlJ l-0 + 2(l-0)2 ^ + 2(l-0)3 \a l3 ) ' 

where £ = tt {n\o\ + n 2 cr|) /6. Comparison with computer simulations for binary molecular 
hard sphere mixtures have shown that the Carnahan-Starling equation ( |İ~5|) is accurate in 
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most of the fluid region, although it fails for high densities and for large diameter ratios 
fT7|| . Given the values considered in our simulations, we expect t hat the approximation ( |Î5|) 
turns out to be quite accurate to evaluate the pair correlation function Xij- in terms of the 
reduced distributions /*, the Enskog kinetic equations become 



<1A 

2 <9v 



where (* = Çi/nvoaf 2 an d J*j = (^o/ nn î (T i2) Jij are given, respectively, by 



6 3=1 J 



(16) 



(17) 



XjXij Vır) I dw * 2 1 da e ^ ' g * u ^ a ' g * 2 ^ 



(18) 



where g^ 2 = gı^Ato) \ — i v o/ v oi) 2 — T/(Tifiji), with j ^ z, is the square of the ratio of the 



thermal velocity to that for species i, and ı>oi = y2Tj/r?v 

In the dimensionless form ( |Î6"D the Enskog equations are time independent. The pair of 
coupled equations C |X6]) must be solved self-consistently with the expressions for the cooling 
rates in Eq. (||) to determine /* and Q = Q — (*■ The temperature T (t) is then obtained 
from the known cooling rate by solving the first of Eqs. (§), and the distribution functions 
f- 1 ' (v, t) are fully determined. The kinetic temperature for each species is obtained from 
Eq. as 

2 



T l (t)=T(t)-^ Jdv*v* 2 f*(v*). 



(19) 



As anticipated, Ti(t) oc T (t) in the HCS and 7 = Tı/T 2 becomes 



7 



(20) 



/i 2 ı J dv*v* 2 fö (v*)' 

This is essentially the approach used in numerical Monte Carlo solution || to the Enskog 
equation. 
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In practice, only approximate solutions for the HCS are possible (an exception is a recent 
exact result for a one dimensional Maxwell model ||18|| ) and a different approach is followed. 
First, the solution is represented as a series in velocity polynomials, with the leading terms 
given by 



/ \ \ 3/2 
f*( v *) "M 7) ^ 



1 + | (A 2 ,* 4 - 5A.,* 2 + Ş) 



(21) 



Thus, the weight function (Gaussian) for each species is chosen to be scaled relative to the 
thermal velocity for that species, introducing explicitly the unknown species temperatures. 
The coefficients c% measure the deviation of /* from the chosen reference Gaussians. The 
cooling rates are now calculated as explicit functions of A, and q from Eq. (|Î7İ). With these 
known, the Enskog equations can be solved to determine q as functions of Aj by substitution 
of (pip into the Enskog equations, taking the v 4 moment of those equations, and retaining 
terms up through linear in q. Finally, the Aj are determined from the consistency condition 
for the HCS, £* = Q. The detailed results for q and Aj as functions of the fluid parameters 
are given in Refs. H and H and will not be repeated here. 



IV. COMPARISON OF THEORY AND SIMULATION 



The approximation (pî|) provides detailed predictions for the species temperatures as 
functions of the mass ratio, size ratio, composition, density, and restitution coefficients. The 
quality of this approximate solution to the Enskog equations has been recently confirmed 
by direct Monte Carlo simulation of those equations över a wide range of the parameter 
space ||. Specifically, the parameter space över which the solution (21) has been verified 



is the mass ratio mı/ma, the concentration ratio nı/712, the ratio of diameters <Jı/(J2, the 
reduced density ncrf 2 , and the (common) restitution coefficient a = an = 0:22 = 0:1.2. 
However, uncertainties remain regarding the accuracy of the Enskog equations themselves. 
An appropriate means to study the concept of the HCS and associated different partial 
temperatures, as well as to study the domain of validity of the Enskog kinetic theory is via 
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MD simulations. Since the parameter space here is quite large the tests of the theory and 
concepts are quite stringent. 

Two different values of the solid volume fractions (f) have been considered here, <fi = 0.1 
and 4> = 0.2, both representing a moderately dense fluid. Ali coefficients of restitution were 
set equal and two values considered, a = 0.8 and a = 0.95, both representing moderately 
strong dissipation. The temperature ratio 7 in the HCS has been studied for three cases in 
each state. In the first case (case I) 7 is determined as a function of the mass ratio 777.1/777.2 
for ö\jo2 = 4>ı/4>2 — 1- The second case (case II) determines 7 as a function of size ratio 
o~ı/o~2 for mı/1712 = (pı/4>2 = 1, while the third case (case III) determines 7 as a function of 
composition <f>\/4>2 for mı/rri2 = 8 and a\jaı — 2. 

The granular system under consideration does not contain external force fields, and thus 
the particles travel in straight-line trajectories between collisions. Correspondingly, an event- 
driven algorithm is employed in the MD simulations. The simulated particles are modeled 
as inelastic, frictionless hard spheres (i. e. collisions are both binary and instantaneous) 
moving in a three-dimensional space with standard periodic boundaries. The initial partide 
velocities are uniformly distributed about a zero mean, regardless of partide size. These 
velocities are then adjusted to ensure that the net system momentum is zero. 

As indicated in Fig. [I], the system reaches a steady value for the temperature ratio 
7 within 10 collisions per partide for a wide class of initial conditions. However, it is 
known that the HCS is unstable to long wavelength perturbations so that spontaneous 
deviations from the HCS occur at long times. To assure that 7 is measured in the HCS, the 
temperature T (t) is monitored as a function of time to determine if the predicted cooling 
from the scaling form ([§]) (Haff's law) is maintained [jî9|1 . In order to keep the computational 
time reasonable for each of the simulations (about 1 hour), the total number of particles 
was kept constant at İV = 1000 for ali simulations. Data from the first 10 collisions per 
partide (or 10000 total collisions) was used to determine the slope of T (t). Specifically, 
the ln(T(t)/T(0)) was sampled 1000 times during this initial portion of the simulation. 
Somewhat surprisingly, a smooth linear decrease in ln(T(£)/T(0)) was observed throughout 
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the 1000 samples, and a linear regression analysis was employed to evaluate the slope of 
the Haff's law plot. Following evaluation of the slope of ln(T(t)/T(0)), collection of the 
7 = Tı(t)/T2(t) data commenced. This data collection period involved as many as 200 
additional collisions per partide (or 200000 total collisions) including 50000 equally spaced 
measurements of 7. The phrase "as many as" refers to the fact that the data collection 
would cease (with fewer than 50000 measurements of the energy ratio) if the measured 
value of ln(T(i)/T(0)) deviated from the expected value of ln(T(t)/T(0)) by more than 5%. 
Violation of the Haff's law restriction occurred frequently when the mass ratio 777.1/777.2 was 
greater than 4. Additionally, simulations of equal mass particles (mı/7772 = 1) violated the 
Haff's law restriction when a = 0.8 and = 0.2. 

Figüre 0shows the results for case I, 7 as a function of mass ratio. The symbols represent 
the simulation data where the circles are for a = 0.95 and the triangles are for a = 0.8. 
In addition, open (solid) symbols correspond to = 0.1 (0 = 0.2). The simulation values 
reported represent the average from three identical simulations, with a standard deviation 
typically less than 3%. The Enskog prediction of the previous section is given by the solid 
lines (the theory does not predict any dependence on in this case). The agreement between 
the theory and simulation is seen to be quite good at a = 0.95, över the whole range of mass 
ratios. The agreement is also quite good at a = 0.8 and = 0.1. However, systemmatic 
deviations from the Enskog theory for large mass ratios are obtained in the simulations at 
= 0.2. 

Figüre ^| shows the results for case II, 7 as a function of size ratio. The notation is the 
same as in Fig. ^| where now the solid line refers to = 0.1 while the dashed line is for 
= 0.2. The agreement for both a = 0.95 and a = 0.8 is quite good at = 0.1, except for 
the largest size ratio at a = 0.8. The density dependence of the theory is weaker than that 
from the simulation, and large differences are observed at = 0.2. 

Figüre ^ shows the results for case III, 7 as a function of composition. We observe 
that both the theory and simulation predicts a very weak inffuence of composition on the 
temperature ratio. In addition, the trends are similar to those of Figs. || and |[ Good 
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agreement is obtained for a = 0.95 at both = 0.1 and <fi = 0.2. At stronger dissipation 
there is a strong density dependence in the simulation that is not reproduced by the theory. 



V. DRIVEN SYSTEMS AND EXPERIMENTS 

The existence and details of different temperatures for each species in a HCS is now well 
established by kinetic theory and simulation. Related experiments |ÎÖ1 , |Î1İ and simulations 



İ4İÎ3H on driven steady states also show different temperatures, but the detailed dependence 
on the control parameters appears to be different. The driven steady states are achieved 
from external forces that do work at the same rate as collisional cooling. In the experiments 
this is accomplished by vibrating the system so that it is locally driven at the walls. Far 
from these walls a steady state is studied whose properties are presumed to be insensitive 
to the details of the driving forces. The velocities of the particles can be measured using 
high speed photography [T(J or positron emission partide tracking [DJ. The objective of 



this section is to explore similarities and differences between the temperature ratios for a 
binary mixture in the HCS and in a driven steady state. 

As a first analysis, a homogeneously driven steady state is considered. This does not 
correspond directly to any experimental driving source but has been considered extensively 
as a representation of driven systems for the one component fluid ||20|| . In this case a 
uniform external non-conservative force, frequently referred to as a "thermostat" , is applied 
to compensate for collisional cooling. Two types of thermostats are considered here. One 
is a deterministle Gaussian thermostat widely used in nonequilibrium molecular dynamics 
simulation for molecular fluids The force is similar to a Stokes law drag force, linear 



in the velocity, but with the opposite sign so that it heats rather than cools the system. 
The "friction" constant can be chosen to exactly compensate for collisional cooling. At the 
level of kinetic theory, the introduetion of such an external force leads to a steady state 
equation which is identical to Eq. (fTüp. It is easily confirmed that the same is true at the 
level of the Liouville equation in appropriate dimensionless variables. Thus, there is an exact 
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correspondence between the HCS and this type of driven steady state and, in particular, the 
dependence of 7 on the control parameters is the same. 

A second method of driving the system homogeneously is by means of a stochastic 
Langevin force representing Gaussian white noise . This force for each species is written 
as Ti = rrii£i, where the covariance of the stochastic acceleration is 

(&a(0&?(O> = ZDÖijÖapÖ (t - t') (22) 

Note that the covariance for the random accelerations is taken to be the same for each species 



PJ23||. This force induces a diffusion in velocity space, with diffusion coefficient D. At the 
level of kinetic theory this leads to an additional source represented by a Fokker-Planck 
collision operatör in addition to the Enskog collision operatör. The steady state Enskog 



equations then take the form 



d^ 2 



= £J« vl/JVf +D[^) if. (23) 



Multiplying by rriivj/2 and integrating gives the relationship of D to the cooling rates d, 
i. e., D = ( i T i /2mi. This in turn implies the steady state condition 

T To 

Cı— = C2— (24) 

The cooling rates are no longer equal, as for the HCS, and the dependence of the tempera- 
tures on the control parameters will be different as well. 

The procedure for determining the temperatures for the stochastically driven steady 
state is the same as that described in Sec. |TT[ The steady state distribution is represented 
as an expansion of the form fl2~l~|) and the coefficients are now determined from moments 
of the set (p3|). The cooling rates are then determined from this solution using Eq. ([TTD , 
and the condition ( p4]) gives an equation for the temperature ratio 7. Figures ^ |^, and |^ 
illustrate the differences between the HCS and the stochastic steady state for a = 0.6 and 
0.8. The solid lines are the results for the HCS while the dashed lines are the results for the 
driven steady state. The dependence of 7 on mass ratio is shown in Fig. || for = and 
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0ı/02 = crı/(T2 = 1. This dependence is seen to be considerably stronger in the driven state. 
The dependence on composition is shown in Fig. |6| for = 0, 1711/1712 = 2 and a\jaı = 1. 
Finally the dependence on overall packing fraction is shown in Fig. ^ for 0ı/02 = 1 and 
mı/1712 = o\jo2 = 2. In this last case the effect of increased density is greater for the HCS 
than for the driven steady state. 

The HCS and homogeneously driven steady states are seen to be qualitatively similar, 
with only quantitative differences. It remains to understand their relationship to locally 
driven wall forces. An example is described for the Boltzmann equation in the Appendbc. 
There, the boundary condition is a sawtooth vibration of one wall such that every partide 
encountering the wall has a reflected speed increased by twice the velocity of the wall in the 
component normal to the wall. The steady state condition is considerably more complex than 
for the HCS or the homogeneously driven steady states. In the limit that the wall velocity is 
large compared to the thermal velocities of each species the condition (p4j) is recovered. This 
suggests that the results obtained from this condition are plausible first approximations for 
qualitative comparisons with experimental results |3|. However, the detailed nature of the 
driven state requires further characterization before quantitative conclusions can be drawn. 
This is suggested by the study of a driven state in the absence of gravity [2~3] where the 
system is found to be well-described by hydrodynamics away from the wall, but the steady 
state is strongly inhomogeneous. 



VI. DISCUSSION 

The primary results of this study are two-fold. First, the MD simulations confirm the 
rapid approach to a HCS with two kinetic temperatures determined by a common cooling 
rate. This occurs över a wide range of densities, composition, mass and size ratios, for both 
moderate and strong dissipation. The second result is confirmation of the Enskog kinetic 
theory to provide a quantitative description of this phenomenon for the lower densities 
and weaker dissipation cases. This includes densities well outside the Boltzmann limit and 
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applies throughout the parameter space of mechanical properties. The analysis here is a 
test of the Enskog prediction for the cooling rates, which are essentially transport properties 
(collision rates). The good agreement obtained is further testimony to the utility of this 
remarkable equation for flııids with elastic and inelastic collisions, including mixtures. 

The failure of the Enskog theory at high densities is expected from experience with normal 
fiuids. This is due to multi-particle collisions, including recollision events (ring collisions). 
The latter are expected to be stronger for fiuids with inelastic collisions where colliding pairs 
tend to become more focused. It appears that the range of densities for which the Enskog 
description applies decreases with increasing dissipation. This is the case observed here and 
also elsewhere for the self-diffusion coefhcient ||. The specific mechanism responsible for 
the large discrepancies at high densities and its quantitative prediction remains an open 
problem. 

The magnitude of the difference between the two kinetic temperatures generally increases 
as the mechanical differences increase, although the dependence on volume fraction is weak. 
Also, there is a significant dependence on the inelasticity and total volume fraction. The 
experiments in |]TÜ[ show a similar strong dependence on mass ratio, but no significant 
dependence on inelasticity, total density, or composition. The detailed correspondence be- 
tween the simple model homogeneous states considered here and the locally driven states of 
experiments needs refinement, although generally the same trends are observed ||. 

The hydrodynamics for binary mbrtures of inelastic hard spheres has been derived re- 
cently, including the effects of two kinetic temperatures ||25|| . Although only the overall 
temperature associated with both species serves as a hydrodynamic field, the transport 
coefficients depend on the temperature ratio 7. Since the latter is a function of the composi- 
tion and density, there are additional contributions to the transport coefficients. Differences 
as large as 50% are found for some coefficients. 
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APPENDIX A: LOCAL BOUNDARY CONDITIONS 



The Enskog equations with boundary conditions can be written as [go| 



(d t + vı • V - T Wİ ) /W(r, v ı; t) = J2 Ja [yı\f\t), tf\t)\ , (Al) 

j 

Here, Twi describes interactions of particles of type i with the boundaries 

T W iîP (r, v l5 1) = j s dsö (r - s) || dv'P( Vl , V) |n • v'| /f (r, v'; t) 

-GC-n-vOln-v^^^V!^)} (A2) 

This is similar to the usual Boltzmann collision operatör with the first term representing a 
gain in the population of particles with velocity Vı due to collisions with the wall and the 
second term representing a corresponding loss. The probability density for a velocity vı 
after the collision with the wall, given an incident velocity v' has the form 

P(vı, v') = 0(n • vı)üf(vı, v')e(-n • v). (A3) 

The two theta functions characterize incident particles as directed toward the wall and 
reflected partide directed away from the wall, where the normal n is directed toward the 
interior of the system. The kernel K(v', Vı) characterizes the change in the half space 
velocity distributions at the surface (i. e., outgoing distribution is a linear functional of the 
incoming distribution). Partide number conservation requires 

J dvxP(v h v) = e(-n • v) (A4) 

As an illustration, the form of K(y,Vı) for elastic specular collisions with a wall at rest is 

17 



K s (v 1 ,v') = 5(v 1 -v' + 2(n-v')n) (A5) 

The surface domain S is time dependent, representing vertical oscillations. This can be 
handled using a sawtooth form for the driving, such that every particle encounters the wall 
moving into the system at velocity v w = v w n. Assume the amplitude of the vibration is 
small so that the displacement of the wall can be neglected. Then, it is reasonable to choose 
specular collisions in the local Galillean frame for which the wall is at rest: 

P(vı, v') = 0(n • (vı- Vt0 ))ür(vı, v')0(-n • (v' - v w )), (A6) 

K s {y x y) = 5 (vı - v' + 2 (n • (v' - v w )) n) . (A7) 

The wall collision operatör becomes 

T Wi fW (r, vı; t) = jf dsö (r - s) {©(n • (vı - v w )) J dv'5 (vı — v' + 2n • (v' - v w ) n) 
x0( _n . ( V ' - v w )) |n • (v' - v w )| /«(r, v'; t) 

-0(-n • ( Vl - vj) |n • ( Vl - v w )\ /f } (r, Vl ; t)} (A8) 

The change in the kinetic energy due to the boundary conditions is 

J dr J dv 1 ^rn i vfT Wi fj: 1 \T,\ 1 -,t)= J^ds J dvı [e- (uı) - e^uı)] |n • (vı - v w )| 

x0 (-n-(v 1 -v w ))/f ) (s,v 1 ;t), (A9) 

where ei(vı) and e -(fi) are the kinetic energies for particles coming into and leaving the 
surface. They are given by 

e i( v ı) = ^ m i v ı, ( A1 °) 

ej(uı) = f rfv0(n • (v - v w ))^m iV 2 ö (v - Vl + 2n • ( Vl - vj n) (Ali) 
= ^m, (vı - 2 (n • vı - îi) 2 0(-n • (vı - v M )). (A12) 

Thus, one gets 
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eiivı) - ei(vı) = 2 mi v 2 w ( 1 - ^-^ J (A13) 

Here, is the speed of the wall, taken along the normal n. The rate of energy change due 
to the wall becomes 

2 



n • V! 



x | = ds0(-fi-(v 1 -vj)// 1) (s,v 1 ;t) (A14) 

If the wall speed is much larger than other characteristic velocities, 1 — n • v/v w — > 1, and 
so Eq. (|A14|) becomes approximately 

J dr J dvı-miv{T Wi fP (r,Yı,t) -> m^rtiA, (A15) 

where A is the area of the wall and rii is the density of incident particles of species i. The 
rate of change of the temperature is then 

= —mtvl - OT, (A16) 

dt W w s v ; 



which gives the steady state condition (23) 



Another estimate is given by representing the incident distribution as a Gaussian 
27^m^— / du* 1 - — v* e - ^ = 2n i Am i v w vl i 

V7TJ-00 V U,ü / 

1 [« + ^*) 2 + 2İ (erf « - ^) + 1)) (A17) 



Here u* = v w /vu, v\ = 2Tîj/mj characterizes the temperature of the incident particles and 
t> 2 = characterizes their mean speed toward the wall. For large v* m + v\ the condition 
(|A16 ) is recovered, while for small v* w + v\ it becomes 

W = W VwTli ~ ClTl (A18) 

If Tu m Ti then the steady state condition for the HCS, constant cooling rate, is recovered. 
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FIGURES 

FIG. 1. Time evolution of 7 (t) = Tı{t)/T 2 {t) for = 0.1, aı/a 2 = fa/fa = 1, mı/m 2 = 8 and 
two values of a: a = 0.8 and a = 1. 

FIG. 2. Plot of the temperature ratio T\jT-ı as a function of the mass ratio mı/m2 for 
C1/C2 = 0ı/ 02 = 1, and two different values of a: a = 0.95 (solid line and circles) and a = 0.8 
(solid line and triangles). The lines are the Enskog predictions and the symbols refer to the MD 
simulation results. The open (solid) symbols correspond to = 0.1 (0 = 0.2). 

FIG 3. Plot of the temperature ratio Î1/T2 as a function of the size ratio o\joı for 
7711/777,2 = 0ı/0 2 = I; and two different values of a: a = 0.95 (lines and circles) and a = 0.8 
(lines and triangles). The lines are the Enskog predictions and the symbols refer to the MD sim- 
ulation results. The solid (dashed) lines correspond to = 0.1 (0 = 0.2) while the open (solid) 
symbols correspond to = 0.1 (0 = 0.2). 

FIG 4. Plot of the temperature ratio Î1/T2 as a function of composition 0ı/0 2 for 7771/7772 = 8, 
cı/02 = 2, and two different values of a: a = 0.95 (lines and circles) and a = 0.8 (lines and 
triangles). The lines are the Enskog predictions and the symbols refer to the MD simulation 
results. The solid (dashed) lines correspond to = 0.1 (0 = 0.2) while the open (solid) symbols 
correspond to = 0.1 (0 = 0.2). 

FIG 5. Plot of the temperature ratio Tı/T 2 as a function of the mass ratio 7771/7772 for = 0, 
C1/C2 = 0ı/02 = 1, and two different values of a: a = 0.8 and a = 0.6. The solid lines are the 
results for the HCS while the dashed lines are the results for the driven steady state achieved from 
the stochastic thermostat. 

FIG 6. Plot of the temperature ratio Tı/T 2 as a function of composition 0ı/0 2 for = 0, 
7771/777,2 = 2, o\j 02 = 1, and two different values of a: a = 0.8 and a = 0.6. The solid lines are 
the results for the HCS while the dashed lines are the results for the driven steady state achieved 
from the stochastic thermostat. 
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FIG. 7. Plot of the relative temperature ratio 7(0, (f))/^f(a, 0) as a function of the total solid 
volume fraction (f> for m\/m2 = 01/02 = 2, fa/fa = 1, and two different values of a: a = 0.8 and 
a = 0.6. The solid lines are the results for the HCS while the dashed lines are the results for the 
driven steady state achieved from the stochastic thermostat. 
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